####################################################################################
### Political leadership has limited impact on fossil fuel taxes and subsidies   ###
### Figures and tables in the supporting information                             ###
### This version: October 29, 2022                                               ###
####################################################################################

#---------------------------------------------------#
# Preliminaries                                     #
#---------------------------------------------------#

rm(list = ls())

#List of packages for session
.packages = c("tidyverse","readstata13", "gridExtra", "lubridate", "estimatr", "texreg",
              "survival", "survminer", "patchwork", "ggplot2")

#Install CRAN packages (if not already installed)
.inst = .packages %in% installed.packages()
if(length(.packages[!.inst]) > 0) install.packages(.packages[!.inst])

#Loading packages into session 
lapply(.packages, require, character.only=TRUE)

# Set Working Directory to wherever files are downloaded
#setwd("")

#---------------------------------------------------#
# Loading Data                                      #
#---------------------------------------------------#

#Main datasets for regression analysis
load("ffs-leaders-2022-10-29.RData")

#Output from the simulations of null distributions
load("complete_null_distributions_oct22.RData")

#Dataset with changes and covariates
leaders_taxchange <- read.csv("leaders_change_newtenure.csv")

#---------------------------------------------------#
# Preparing Data                                    #
#---------------------------------------------------#

#Leader individual covariates

leadership_covariates <- newdat %>%
  distinct(obsid, .keep_all = T) %>%
  dplyr::select(obsid, leader, leader_id, country2, gender, leader_age, unieducation, execrlc)

#Dataset with changes and covariates
leaders_taxchange <- leaders_taxchange %>%
  dplyr::select(obsid, leader, country2, tenure_months, tenure_months_effective,
                mean_tax_entry, mean_tax_exit, diff_tax, rate_change_tax_eff, 
                presidential_democracy, fixedprice, persistent, oil, democracy, 
                pct_crisis_tenure, avg_inflation)

leaders_taxchange_cov <- leaders_taxchange %>%
  left_join(leadership_covariates, by = c("obsid","leader", "country2")) %>%
  mutate(leader_age = case_when(leader_id == "Nouri Abusahmain-LBY" ~ 58,
                                TRUE ~ as.numeric(leader_age)),
         leader_age_cat = case_when(leader_age <= 50 ~ "50 and Younger",
                                    leader_age > 50 & leader_age <= 60 ~ "51 to 60",
                                    leader_age > 60 ~ "61 and Older"),
         democracy = case_when(democracy >= 0.5 ~ "Democracy",
                               TRUE ~ "Autocracy"),
         presidential_democracy = case_when(democracy == "Democracy" & presidential_democracy >= 0.5 ~ "Presidential",
                                            democracy == "Democracy" & presidential_democracy < 0.5 ~ "Parliamentary",
                                            TRUE ~ "Autocracy"),
         unieducation = case_when(unieducation == 1 ~ "College-Level Education",
                                  unieducation == 0 ~ "Non-College-Level Education"),
         ideological_orientation = case_when(execrlc == 1 ~ "Right",
                                             execrlc == 2 ~ "Center",
                                             execrlc == 3 ~ "Left",
                                             TRUE ~ "Non-Applicable"),
         persistent_subsidizer = case_when(persistent == 1 ~ "Persistent Subsidizer",
                                           TRUE ~ "Persistent Taxer"),
         coe = ifelse(leader %in% c("Narendra Modi", "Michelle Bachelet", "Paul Kagame",
                                    "Hasina Wazed", "Bambang Yudhoyono", "Elbegdorj", "Calderon",
                                    "Bharrat Jagdeo", "Helen Clark"), 1, 0),
         oil_exporter = case_when(oil > 0 ~ "Oil Exporter",
                                  TRUE ~ "Non-Oil Importer"),
         pct_crisis_cat = if_else(pct_crisis_tenure >= 0.5, "Persistent Crisis", "Non-Persistent Crisis"))

####################################################
######## Supporting Information -- Figure 1 ########                     
####################################################

## Figure S1. Kernel Density Plots by Subgroup. Democracies versus Autocracies

bmgapchange_regime_plot <- leaders_taxchange_cov %>%
  filter(!is.na(democracy)) %>%
  ggplot(aes(x = rate_change_tax_eff, fill = democracy)) + 
  geom_density(alpha = 0.2) +
  theme_classic() +
  xlim(-0.08, 0.08) +
  xlab("Average Monthly Change in Fuel Tax (USD)") +
  ylab("Count") +
  scale_fill_manual(values = c("#2b7cb3","#ff7f50"), name='Regime Type') +
  theme(legend.key.size = unit(0.3, 'cm'))

bmgapchange_regime_plot

####################################################
######## Supporting Information -- Figure 2 ########                     
####################################################

## Figure S2. Kernel Density Plots by Subgroup. Parliamentary versus Presidential

bmgapchange_presidential_plot <- leaders_taxchange_cov %>%
  filter(!is.na(presidential_democracy)) %>%
  filter(presidential_democracy != "Autocracy") %>%
  ggplot(aes(x = rate_change_tax_eff, fill = presidential_democracy)) + 
  geom_density(alpha = 0.2) +
  theme_classic() +
  xlim(-0.08, 0.08) +
  xlab("Average Monthly Change in Fuel Tax (USD)") +
  ylab("Count") +
  scale_fill_manual(values = c("#2b7cb3","#ff7f50", "#228b22"), name='Democracy Type')  +
  theme(legend.key.size = unit(0.3, 'cm'))

bmgapchange_presidential_plot

####################################################
######## Supporting Information -- Figure 3 ########                     
####################################################

## Figure S3. Kernel Density Plots by Subgroup. Oil Exporters versus Oil Importers

bmgapchange_exporter_plot <- leaders_taxchange_cov %>%
  filter(!is.na(oil_exporter)) %>%
  ggplot(aes(x = rate_change_tax_eff, fill = oil_exporter)) + 
  geom_density(alpha = 0.2) +
  theme_classic() +
  xlim(-0.1, 0.1) +
  xlab("Average Monthly Change in Fuel Tax (USD)") +
  ylab("Count") +
  scale_fill_manual(values = c("#2b7cb3","#ff7f50"), name='Oil Exporters Status')  +
  theme(legend.key.size = unit(0.3, 'cm'))

bmgapchange_exporter_plot

####################################################
######## Supporting Information -- Figure 4 ########                     
####################################################

## Figure S4. Kernel Density Plots by Subgroup. Subsidizers versus Non Subsidizers

bmgapchange_subsidizer_plot <- leaders_taxchange_cov %>%
  filter(!is.na(persistent)) %>%
  ggplot(aes(x = rate_change_tax_eff, fill = persistent_subsidizer)) + 
  geom_density(alpha = 0.2) +
  theme_classic() +
  xlim(-0.1, 0.1) +
  ylim(0,150) +
  xlab("Average Monthly Change in Fuel Tax (USD)") +
  ylab("Count") +
  scale_fill_manual(values = c("#2b7cb3","#ff7f50"), name='Subsidizer Status')  +
  theme(legend.key.size = unit(0.3, 'cm'))

bmgapchange_subsidizer_plot

####################################################
######## Supporting Information -- Figure 5 ########                     
####################################################

## Figure S5. Kernel Density Plots by Subgroup. Crisis versus Non-Crisis

bmgapchange_crisis_plot <- leaders_taxchange_cov %>%
  filter(!is.na(pct_crisis_cat)) %>%
  ggplot(aes(x = rate_change_tax_eff, fill = pct_crisis_cat)) + 
  geom_density(alpha = 0.2) +
  theme_classic() +
  xlim(-0.08, 0.08) +
  xlab("Average Monthly Change in Fuel Tax (USD)") +
  ylab("Count") +
  scale_fill_manual(values = c("#2b7cb3","#ff7f50"), name='Persistent Crisis')  +
  theme(legend.key.size = unit(0.3, 'cm'))

bmgapchange_crisis_plot

####################################################
######## Supporting Information -- Figure 6 ########                     
####################################################

## Figure S6. Mean monthly change in net implicit taxes for leaders

bmgapchange_age_plot <- leaders_taxchange_cov %>%
  filter(!is.na(leader_age_cat)) %>%
  ggplot(aes(x = rate_change_tax_eff, fill = leader_age_cat)) + 
  geom_density(alpha = 0.2) +
  theme_classic() +
  xlim(-0.08, 0.08) +
  xlab("Average Monthly Change in Fuel Tax (USD)") +
  ylab("Count") +
  scale_fill_manual(values = c("#2b7cb3","#ff7f50", "#7ee6a3"), name='Age Cohort') +
  theme(legend.key.size = unit(0.3, 'cm'))

bmgapchange_age_plot

#Gender

bmgapchange_gender_plot <- leaders_taxchange_cov %>%
  filter(!is.na(gender)) %>%
  ggplot(aes(x = rate_change_tax_eff, fill = gender)) + 
  geom_density(alpha = 0.2) +
  theme_classic() +
  xlim(-0.08, 0.08) +
  xlab("Average Monthly Change in Fuel Tax (USD)") +
  ylab("Count") +
  scale_fill_manual(values = c("#2b7cb3","#ff7f50"), name='Gender') +
  theme(legend.key.size = unit(0.3, 'cm'))

bmgapchange_gender_plot

#College Education

bmgapchange_unieducation_plot <- leaders_taxchange_cov %>%
  filter(!is.na(unieducation)) %>%
  ggplot(aes(x = rate_change_tax_eff, fill = unieducation)) + 
  geom_density(alpha = 0.2) +
  theme_classic() +
  xlim(-0.08, 0.08) +
  xlab("Average Monthly Change in Fuel Tax (USD)") +
  ylab("Count") +
  scale_fill_manual(values = c("#2b7cb3","#ff7f50"), name='College Education') +
  theme(legend.key.size = unit(0.3, 'cm'))

bmgapchange_unieducation_plot

#Ideological Orientation

bmgapchangeideological_orientation_plot <- leaders_taxchange_cov %>%
  filter(!is.na(ideological_orientation)) %>%
  ggplot(aes(x = rate_change_tax_eff, fill = ideological_orientation)) + 
  geom_density(alpha = 0.2) +
  theme_classic() +
  xlim(-0.08, 0.08) +
  xlab("Average Monthly Change in Fuel Tax (USD)") +
  ylab("Count") +
  scale_fill_manual(values = c("#2b7cb3","#ff7f50", "#7ee6a3",
                               "#FFEE58"), name='Ideological Orientation') +
  theme(legend.key.size = unit(0.3, 'cm'))

bmgapchangeideological_orientation_plot

####################################################
######## Supporting Information -- Figure 7 ########                     
####################################################

## Figure S7. Fossil Fuel Subsidies Reforms Survival Plots

#-----------------------------------------------------------------------------#
#    Reforms for PS countries using Nominal LCU price                         #
#-----------------------------------------------------------------------------#

# Monthly change in Nominal LCU price for persistent subsidizers
monthly_ps <- censdat %>%
  filter(persistent == 1) %>%
  arrange(country2, year_month) %>%
  group_by(country2) %>%
  mutate(price.1 = lag(price, n = 1),
         usd.1 = lag(price_usd_2015, n = 1),
         bmg.1 = lag(bmgap2015adj, n = 1),
         delta.nom = price - price.1,
         d.pct.nom = (price/price.1) - 1,
         delta.bmg = bmgap2015adj - bmg.1,
         delta.usd = price_usd_2015 - usd.1)

# Coding reform as any nominal change above 10% threshold
monthly_ps$reform <- 1*(monthly_ps$delta.nom > 0 & monthly_ps$d.pct.nom > 0.1)
monthly_ps$reform[is.na(monthly_ps$reform)] <- 0

# Creating reform-only list
monthly_ps_reforms <- list()
nom_reforms <- list()
usd_reforms <- list()
bmg_reforms <- list()

monthly_ps_select <- monthly_ps[,c("country2","year_month","date","leader_id","price",
                                   "price_usd_2015","bmgap2015adj", "price.1","usd.1",
                                   "bmg.1","d.pct.nom","reform")]

for(i in 1:sum(monthly_ps$reform)){
  # for(i in 1:2){
  
  # Identifying reforms and storing each reform as its own dataframe
  monthly_ps_reforms[[i]] <- monthly_ps_select[monthly_ps_select$reform==1,][i,]
  
  # Binding each reform dataframe with all months that follow
  monthly_ps_reforms[[i]] <- rbind(monthly_ps_reforms[[i]], monthly_ps_select[monthly_ps_select$country2==monthly_ps_reforms[[i]]$country2 & monthly_ps_select$date > monthly_ps_reforms[[i]]$date,])
  
  # Reform date
  monthly_ps_reforms[[i]]$reform_start <- monthly_ps_reforms[[i]]$date[1] 
  
  # Pre-reform prices: nominal, real USD, BMG
  monthly_ps_reforms[[i]]$pre.reform.price <- monthly_ps_reforms[[i]]$price.1[1]
  monthly_ps_reforms[[i]]$pre.reform.usd <- monthly_ps_reforms[[i]]$usd.1[1]
  monthly_ps_reforms[[i]]$pre.reform.bmg <- monthly_ps_reforms[[i]]$bmg.1[1]
  
  # Failure for each type
  monthly_ps_reforms[[i]]$reform_nom_fail <- 1*(monthly_ps_reforms[[i]]$price <= monthly_ps_reforms[[i]]$pre.reform.price)
  monthly_ps_reforms[[i]]$reform_usd_fail <- 1*(monthly_ps_reforms[[i]]$price_usd_2015 <= monthly_ps_reforms[[i]]$pre.reform.usd)
  monthly_ps_reforms[[i]]$reform_bmg_fail <- 1*(monthly_ps_reforms[[i]]$bmgap2015adj <= monthly_ps_reforms[[i]]$pre.reform.bmg)
  
  # Spell length
  monthly_ps_reforms[[i]]$spell <- as.numeric(difftime(monthly_ps_reforms[[i]]$date, monthly_ps_reforms[[i]]$reform_start, units = "days")) / 365.25
  
  # Reducing to start and end (Nominal)
  nom_reforms[[i]] <- monthly_ps_reforms[[i]]
  
  if(sum(nom_reforms[[i]]$reform_nom_fail) > 0){
    nom_reforms[[i]] <- nom_reforms[[i]][nom_reforms[[i]]$reform_nom_fail==1,][1,]
    nom_reforms[[i]]$status <- 1
  } else {
    nom_reforms[[i]] <- nom_reforms[[i]][1,]
    nom_reforms[[i]]$spell <- as.numeric(difftime(as.Date("2015-06-01"), nom_reforms[[i]]$reform_start, units = "days")) / 365.25
    nom_reforms[[i]]$status <- 0
  }
  
  # Reducing to start and end (USD)
  usd_reforms[[i]] <- monthly_ps_reforms[[i]]
  
  if(sum(usd_reforms[[i]]$reform_usd_fail) > 0){
    usd_reforms[[i]] <- usd_reforms[[i]][usd_reforms[[i]]$reform_usd_fail==1,][1,]
    usd_reforms[[i]]$status <- 1
  } else {
    usd_reforms[[i]] <- usd_reforms[[i]][1,]
    usd_reforms[[i]]$spell <- as.numeric(difftime(as.Date("2015-06-01"), usd_reforms[[i]]$reform_start, units = "days")) / 365.25
    usd_reforms[[i]]$status <- 0
  }
  
  # Reducing to start and end (BMG)
  bmg_reforms[[i]] <- monthly_ps_reforms[[i]]
  
  if(sum(bmg_reforms[[i]]$reform_bmg_fail) > 0){
    bmg_reforms[[i]] <- bmg_reforms[[i]][bmg_reforms[[i]]$reform_bmg_fail==1,][1,]
    bmg_reforms[[i]]$status <- 1
  } else {
    bmg_reforms[[i]] <- bmg_reforms[[i]][1,]
    bmg_reforms[[i]]$spell <- as.numeric(difftime(as.Date("2015-06-01"), bmg_reforms[[i]]$reform_start, units = "days")) / 365.25
    bmg_reforms[[i]]$status <- 0
  }
}

nom_all_reforms <- bind_rows(nom_reforms, .id = "column_label")
usd_all_reforms <- bind_rows(usd_reforms, .id = "column_label")
bmg_all_reforms <- bind_rows(bmg_reforms, .id = "column_label")

#-----------------------------------------------------------------------------#
#    Reforms for all FPS countries using Nominal LCU price                    #
#-----------------------------------------------------------------------------#

# Monthly change in Nominal LCU price for fixed price
monthly_fps <- censdat %>%
  filter(fixedprice == 1) %>%
  arrange(country2, year_month) %>%
  group_by(country2) %>%
  mutate(price.1 = lag(price, n = 1),
         usd.1 = lag(price_usd_2015, n = 1),
         bmg.1 = lag(bmgap2015adj, n = 1),
         delta.nom = price - price.1,
         d.pct.nom = (price/price.1) - 1,
         delta.bmg = bmgap2015adj - bmg.1,
         delta.usd = price_usd_2015 - usd.1)

# Coding reform as any nominal change above 10% threshold
monthly_fps$reform <- 1*(monthly_fps$delta.nom > 0 & monthly_fps$d.pct.nom > 0.1)
monthly_fps$reform[is.na(monthly_fps$reform)] <- 0

# Creating reform-only list
monthly_fps_reforms <- list()
nom_reforms_fps <- list()
usd_reforms_fps <- list()
bmg_reforms_fps <- list()

monthly_fps_select <- monthly_fps[,c("country2","year_month","date","leader_id","price",
                                     "price_usd_2015","bmgap2015adj", "price.1","usd.1",
                                     "bmg.1","d.pct.nom","reform")]

for(i in 1:sum(monthly_fps$reform)){
  # for(i in 1:2){
  
  # Identifying reforms and storing each reform as its own dataframe
  monthly_fps_reforms[[i]] <- monthly_fps_select[monthly_fps_select$reform==1,][i,]
  
  # Binding each reform dataframe with all months that follow
  monthly_fps_reforms[[i]] <- rbind(monthly_fps_reforms[[i]], monthly_fps_select[monthly_fps_select$country2==monthly_fps_reforms[[i]]$country2 & monthly_fps_select$date > monthly_fps_reforms[[i]]$date,])
  
  # Reform date
  monthly_fps_reforms[[i]]$reform_start <- monthly_fps_reforms[[i]]$date[1] 
  
  # Pre-reform prices: nominal, real USD, BMG
  monthly_fps_reforms[[i]]$pre.reform.price <- monthly_fps_reforms[[i]]$price.1[1]
  monthly_fps_reforms[[i]]$pre.reform.usd <- monthly_fps_reforms[[i]]$usd.1[1]
  monthly_fps_reforms[[i]]$pre.reform.bmg <- monthly_fps_reforms[[i]]$bmg.1[1]
  
  # Failure for each type
  monthly_fps_reforms[[i]]$reform_nom_fail <- 1*(monthly_fps_reforms[[i]]$price <= monthly_fps_reforms[[i]]$pre.reform.price)
  monthly_fps_reforms[[i]]$reform_usd_fail <- 1*(monthly_fps_reforms[[i]]$price_usd_2015 <= monthly_fps_reforms[[i]]$pre.reform.usd)
  monthly_fps_reforms[[i]]$reform_bmg_fail <- 1*(monthly_fps_reforms[[i]]$bmgap2015adj <= monthly_fps_reforms[[i]]$pre.reform.bmg)
  
  # Spell length
  monthly_fps_reforms[[i]]$spell <- as.numeric(difftime(monthly_fps_reforms[[i]]$date, monthly_fps_reforms[[i]]$reform_start, units = "days")) / 365.25
  
  # Reducing to start and end (Nominal)
  nom_reforms_fps[[i]] <- monthly_fps_reforms[[i]]
  
  if(sum(nom_reforms_fps[[i]]$reform_nom_fail) > 0){
    nom_reforms_fps[[i]] <- nom_reforms_fps[[i]][nom_reforms_fps[[i]]$reform_nom_fail==1,][1,]
    nom_reforms_fps[[i]]$status <- 1
  } else {
    nom_reforms_fps[[i]] <- nom_reforms_fps[[i]][1,]
    nom_reforms_fps[[i]]$spell <- as.numeric(difftime(as.Date("2015-06-01"), nom_reforms_fps[[i]]$reform_start, units = "days")) / 365.25
    nom_reforms_fps[[i]]$status <- 0
  }
  
  # Reducing to start and end (USD)
  usd_reforms_fps[[i]] <- monthly_fps_reforms[[i]]
  
  if(sum(usd_reforms_fps[[i]]$reform_usd_fail) > 0){
    usd_reforms_fps[[i]] <- usd_reforms_fps[[i]][usd_reforms_fps[[i]]$reform_usd_fail==1,][1,]
    usd_reforms_fps[[i]]$status <- 1
  } else {
    usd_reforms_fps[[i]] <- usd_reforms_fps[[i]][1,]
    usd_reforms_fps[[i]]$spell <- as.numeric(difftime(as.Date("2015-06-01"), usd_reforms_fps[[i]]$reform_start, units = "days")) / 365.25
    usd_reforms_fps[[i]]$status <- 0
  }
  
  # Reducing to start and end (BMG)
  bmg_reforms_fps[[i]] <- monthly_fps_reforms[[i]]
  
  if(sum(bmg_reforms_fps[[i]]$reform_bmg_fail) > 0){
    bmg_reforms_fps[[i]] <- bmg_reforms_fps[[i]][bmg_reforms_fps[[i]]$reform_bmg_fail==1,][1,]
    bmg_reforms_fps[[i]]$status <- 1
  } else {
    bmg_reforms_fps[[i]] <- bmg_reforms_fps[[i]][1,]
    bmg_reforms_fps[[i]]$spell <- as.numeric(difftime(as.Date("2015-06-01"), bmg_reforms_fps[[i]]$reform_start, units = "days")) / 365.25
    bmg_reforms_fps[[i]]$status <- 0
  }
}

nom_all_reforms_fps <- bind_rows(nom_reforms_fps, .id = "column_label")
usd_all_reforms_fps <- bind_rows(usd_reforms_fps, .id = "column_label")
bmg_all_reforms_fps <- bind_rows(bmg_reforms_fps, .id = "column_label")

#-----------------------------------------------------------------------------#
#    Reforms for non-PS FPS countries using Nominal LCU price                 #
#-----------------------------------------------------------------------------#

# Monthly change in Nominal LCU price for non-persistent subsidizers fixed price
monthly_npsfps <- censdat %>%
  filter(fixedprice == 1, persistent == 0) %>%
  arrange(country2, year_month) %>%
  group_by(country2) %>%
  mutate(price.1 = lag(price, n = 1),
         usd.1 = lag(price_usd_2015, n = 1),
         bmg.1 = lag(bmgap2015adj, n = 1),
         delta.nom = price - price.1,
         d.pct.nom = (price/price.1) - 1,
         delta.bmg = bmgap2015adj - bmg.1,
         delta.usd = price_usd_2015 - usd.1)

# Coding reform as any nominal change above 10% threshold
monthly_npsfps$reform <- 1*(monthly_npsfps$delta.nom > 0 & monthly_npsfps$d.pct.nom > 0.1)
monthly_npsfps$reform[is.na(monthly_npsfps$reform)] <- 0

# Creating reform-only list
monthly_npsfps_reforms <- list()
nom_reforms_npsfps <- list()
usd_reforms_npsfps <- list()
bmg_reforms_npsfps <- list()

monthly_npsfps_select <- monthly_npsfps[,c("country2","year_month","date","leader_id","price",
                                           "price_usd_2015","bmgap2015adj", "price.1","usd.1",
                                           "bmg.1","d.pct.nom","reform")]

for(i in 1:sum(monthly_npsfps$reform)){
  # for(i in 1:2){
  
  # Identifying reforms and storing each reform as its own dataframe
  monthly_npsfps_reforms[[i]] <- monthly_npsfps_select[monthly_npsfps_select$reform==1,][i,]
  
  # Binding each reform dataframe with all months that follow
  monthly_npsfps_reforms[[i]] <- rbind(monthly_npsfps_reforms[[i]], monthly_npsfps_select[monthly_npsfps_select$country2==monthly_npsfps_reforms[[i]]$country2 & monthly_npsfps_select$date > monthly_npsfps_reforms[[i]]$date,])
  
  # Reform date
  monthly_npsfps_reforms[[i]]$reform_start <- monthly_npsfps_reforms[[i]]$date[1] 
  
  # Pre-reform prices: nominal, real USD, BMG
  monthly_npsfps_reforms[[i]]$pre.reform.price <- monthly_npsfps_reforms[[i]]$price.1[1]
  monthly_npsfps_reforms[[i]]$pre.reform.usd <- monthly_npsfps_reforms[[i]]$usd.1[1]
  monthly_npsfps_reforms[[i]]$pre.reform.bmg <- monthly_npsfps_reforms[[i]]$bmg.1[1]
  
  # Failure for each type
  monthly_npsfps_reforms[[i]]$reform_nom_fail <- 1*(monthly_npsfps_reforms[[i]]$price <= monthly_npsfps_reforms[[i]]$pre.reform.price)
  monthly_npsfps_reforms[[i]]$reform_usd_fail <- 1*(monthly_npsfps_reforms[[i]]$price_usd_2015 <= monthly_npsfps_reforms[[i]]$pre.reform.usd)
  monthly_npsfps_reforms[[i]]$reform_bmg_fail <- 1*(monthly_npsfps_reforms[[i]]$bmgap2015adj <= monthly_npsfps_reforms[[i]]$pre.reform.bmg)
  
  # Spell length
  monthly_npsfps_reforms[[i]]$spell <- as.numeric(difftime(monthly_npsfps_reforms[[i]]$date, monthly_npsfps_reforms[[i]]$reform_start, units = "days")) / 365.25
  
  # Reducing to start and end (Nominal)
  nom_reforms_npsfps[[i]] <- monthly_npsfps_reforms[[i]]
  
  if(sum(nom_reforms_npsfps[[i]]$reform_nom_fail) > 0){
    nom_reforms_npsfps[[i]] <- nom_reforms_npsfps[[i]][nom_reforms_npsfps[[i]]$reform_nom_fail==1,][1,]
    nom_reforms_npsfps[[i]]$status <- 1
  } else {
    nom_reforms_npsfps[[i]] <- nom_reforms_npsfps[[i]][1,]
    nom_reforms_npsfps[[i]]$spell <- as.numeric(difftime(as.Date("2015-06-01"), nom_reforms_npsfps[[i]]$reform_start, units = "days")) / 365.25
    nom_reforms_npsfps[[i]]$status <- 0
  }
  
  # Reducing to start and end (USD)
  usd_reforms_npsfps[[i]] <- monthly_npsfps_reforms[[i]]
  
  if(sum(usd_reforms_npsfps[[i]]$reform_usd_fail) > 0){
    usd_reforms_npsfps[[i]] <- usd_reforms_npsfps[[i]][usd_reforms_npsfps[[i]]$reform_usd_fail==1,][1,]
    usd_reforms_npsfps[[i]]$status <- 1
  } else {
    usd_reforms_npsfps[[i]] <- usd_reforms_npsfps[[i]][1,]
    usd_reforms_npsfps[[i]]$spell <- as.numeric(difftime(as.Date("2015-06-01"), usd_reforms_npsfps[[i]]$reform_start, units = "days")) / 365.25
    usd_reforms_npsfps[[i]]$status <- 0
  }
  
  # Reducing to start and end (BMG)
  bmg_reforms_npsfps[[i]] <- monthly_npsfps_reforms[[i]]
  
  if(sum(bmg_reforms_npsfps[[i]]$reform_bmg_fail) > 0){
    bmg_reforms_npsfps[[i]] <- bmg_reforms_npsfps[[i]][bmg_reforms_npsfps[[i]]$reform_bmg_fail==1,][1,]
    bmg_reforms_npsfps[[i]]$status <- 1
  } else {
    bmg_reforms_npsfps[[i]] <- bmg_reforms_npsfps[[i]][1,]
    bmg_reforms_npsfps[[i]]$spell <- as.numeric(difftime(as.Date("2015-06-01"), bmg_reforms_npsfps[[i]]$reform_start, units = "days")) / 365.25
    bmg_reforms_npsfps[[i]]$status <- 0
  }
}

nom_all_reforms_npsfps <- bind_rows(nom_reforms_npsfps, .id = "column_label")
usd_all_reforms_npsfps <- bind_rows(usd_reforms_npsfps, .id = "column_label")
bmg_all_reforms_npsfps <- bind_rows(bmg_reforms_npsfps, .id = "column_label")


# Drop the large lists
remove(monthly_ps_reforms, monthly_fps_reforms, monthly_npsfps_reforms)

#-----------------------------------------------------------------------------#
# FPS, PS, Non-PS FPS overlap
#-----------------------------------------------------------------------------#

unique(censdat$country2[censdat$persistent==1]) %in% 
  unique(censdat$country2[censdat$fixedprice==1])

unique(censdat$country2[censdat$persistent==0])[which(unique(censdat$country2[censdat$persistent==0]) %in% unique(censdat$country2[censdat$fixedprice==1]))]

#-----------------------------------------------------------------------------#
# In-text mentions to reforms
#-----------------------------------------------------------------------------#

# Annual instances of reform
# "about one reform every 4 country-years"
monthly_ps %>%
  group_by(country2,year) %>%
  summarize(reforms = sum(reform),
            nums = n()) %>%
  nrow() / sum(monthly_ps$reform)

nrow(monthly_ps)

mean(nom_all_reforms$spell)*12
mean(usd_all_reforms$spell)*12
mean(bmg_all_reforms$spell)*12


# Iran example - Khatami 1999

monthly_ps %>%
  filter(country2=="IRN", leader=="Khatami", reform==1, year == 1999) %>%
  dplyr::select(year_month, time_trend, price, price.1, price_usd_2015, usd.1, bmgap2015adj, bmg.1)

usd_all_reforms %>%
  filter(leader_id=="Khatami-IRN") %>% 
  arrange(year_month) %>%
  dplyr::select(year_month, price, price.1, price_usd_2015, pre.reform.usd, reform_usd_fail)

#(as.Date("2002-03-01") - as.Date("1999-04-01"))/30.5

bmg_all_reforms %>%
  filter(leader_id=="Khatami-IRN") %>% 
  arrange(year_month) %>%
  dplyr::select(year_month, price, price.1, bmgap2015adj, pre.reform.bmg, reform_bmg_fail)

#-----------------------------------------------------------------------------#
# Survival plots
#-----------------------------------------------------------------------------#

# PS Reforms (nominal, real USD, BMG)

nom_fit <- survfit(Surv(spell, status) ~ 1, data = nom_all_reforms)
usd_fit <- survfit(Surv(spell, status) ~ 1, data = usd_all_reforms)
bmg_fit <- survfit(Surv(spell, status) ~ 1, data = bmg_all_reforms)

fit <- list(realUSD = usd_fit, nominal = nom_fit, bmg = bmg_fit)

plot_ps <- ggsurvplot(
  fit = fit, combine = TRUE,
  conf.int = FALSE,                    # Add confidence interval
  conf.int.style = "ribbon",            # CI style, use "step" or "ribbon"
  censor = FALSE,                     # Remove censor points
  # risk.table = TRUE,
  # tables.height = 0.2,
  tables.theme = theme_cleantable(),
  title = "Persistent Subsidizers: Duration of large reforms",
  xlab = "Years", 
  ylab = "Overall survival probability", 
  legend.title = "",
  xlim = c(0,10))

# plot_ps$plot <- plot_ps$plot + scale_x_continuous(breaks = seq(0,20,5))
plot_ps$plot

# As three separate plots
plotcols <- ggsci::pal_aaas(palette = c("default"), alpha = 1)(3)

plot3.nom <-ggsurvplot(
  fit = nom_fit, #combine = FALSE,
  color = plotcols[1],
  conf.int = TRUE,                    # Add confidence interval
  conf.int.style = "ribbon",            # CI style, use "step" or "ribbon"
  censor = FALSE,                     # Remove censor points
  risk.table = TRUE,
  tables.height = 0.2,
  tables.theme = theme_cleantable(),
  title = "A: Reform duration, nominal",
  xlab = "Years", 
  ylab = "Overall survival probability", 
  legend = "none",
  risk.table.title = "Number of enduring reforms",
  risk.table.y.text = F,
  xlim = c(0,10)
)

plot3.usd <-ggsurvplot(
  fit = usd_fit, #combine = FALSE,
  color = plotcols[2],
  conf.int = TRUE,                    # Add confidence interval
  conf.int.style = "ribbon",            # CI style, use "step" or "ribbon"
  censor = FALSE,                     # Remove censor points
  risk.table = TRUE,
  tables.height = 0.2,
  tables.theme = theme_cleantable(),
  title = "B: Reform duration, real",
  xlab = "Years", 
  ylab = "", 
  legend = "none",
  risk.table.title = "Number of enduring reforms",
  risk.table.y.text = F,
  xlim = c(0,10)
)

plot3.bmg <- ggsurvplot(
  fit = bmg_fit, #combine = FALSE,
  color = plotcols[3],
  conf.int = TRUE,                    # Add confidence interval
  conf.int.style = "ribbon",            # CI style, use "step" or "ribbon"
  censor = FALSE,                     # Remove censor points
  risk.table = TRUE,
  tables.height = 0.2,
  tables.theme = theme_cleantable(),
  title = "C: Reform duration, price gap",
  xlab = "Years", 
  ylab = "", 
  legend = "none",
  risk.table.title = "Number of enduring reforms",
  risk.table.y.text = F,
  xlim = c(0,10)
)

pl = list(plot3.nom$plot, plot3.usd$plot, plot3.bmg$plot)#, 
#plot3.nom$table, plot3.usd$table, plot3.bmg$table)
margin = theme(plot.margin = unit(c(0.5,1,0.5,1), "cm"))

gridExtra::grid.arrange(grobs = lapply(pl, "+", margin), 
                        ncol=3)

####################################################
######## Supporting Information -- Figure 8 ########                     
####################################################

## Figure S8. Distribution of Government Tenures with Crises

leaders_taxchange_cov_ordered <- leaders_taxchange_cov %>%
  arrange(rate_change_tax_eff) %>% 
  mutate(orders = row_number())

# Color binary
ggplot(leaders_taxchange_cov_ordered %>%
         filter(!is.na(pct_crisis_cat)), 
       aes(x = rate_change_tax_eff, y = orders, 
           col = pct_crisis_cat)) +
  geom_point() +
  labs(x = "Average rate of change in fuel taxes, by leader term", y = "") +
  # scale_colour_distiller(palette = "RdYlGn") +
  scale_color_manual(name = "Persistence indicator: \n>50% of tenure during crisis", values = c("#3d9b43","#d43025")) +
  # scale_alpha_manual(values = c(0.5, 0.5)) +
  # ggrepel::geom_label_repel(data = leaders_taxchange_cov_ordered %>% filter(pct_crisis_cat == "Persistent Crisis"), 
  #                           aes(x = rate_change_tax_eff, y = orders, label = leader)) +
  scale_x_continuous(breaks = seq(-0.05,0.05,0.025)) +
  theme(legend.position = c(0.2,0.85),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        plot.title = element_text(size=22),
        axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=12), 
        axis.line.y = element_blank(),
        axis.text.y = element_blank(), 
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank())

####################################################
######## Supporting Information -- Table S1 ########                     
####################################################

## Creating the list of countries to be included in the sample

## Dropped countries:"AGO" "BLR" "CMR" "COG" "GMB" "GNQ" "KAZ" "LBR" "MMR" "OMN" "RWA" "SDN" "SWZ" "SYR" "TCD" "TJK" "UGA" "UZB" "ZWE"

oneleadercountry <- censdat %>% group_by(country2) %>% summarise(n_lead = n_distinct(obsid)) %>% filter(n_lead > 1)

## Creating new tenure length, leader trend, and country trend variables

censdat <- censdat %>%
  filter(country2 %in% oneleadercountry$country2) %>%
  group_by(country2) %>%
  mutate(time_trend_country = row_number()) %>%
  ungroup() %>%
  group_by(obsid) %>%
  mutate(tenure_months_new = n(),
         time_trend_leader = row_number())

#Table S1. Model results for democracies

#For all models, column 1 has country intercepts and country trends,
#column 2 has country intercepts and leader traits
#column 3 has leader fixed effects

#Table 2: All Democracies

fit.cfe1.dem <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(country2)
                          + as.factor(country2)*time_trend_country, 
                          data=censdat[censdat$democracy==1,],
                          clusters = country2)

fit.lead.traits.dem <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                                 + NY.GDP.MKTP.KD.ZG + percent_gdp
                                 + fuel_income_dependence + VAT.Rate
                                 + leader_age + gender + unieducation + execrlc
                                 + as.factor(country2)
                                 + as.factor(country2)*time_trend_country,
                                 data=censdat[censdat$democracy==1,],
                                 clusters = country2)

fit.lfe1.dem <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(obsid)
                          + as.factor(obsid)*time_trend_leader, 
                          data=censdat[censdat$democracy==1,],
                          clusters = country2)

fit.crisis1.dem <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                             + NY.GDP.MKTP.KD.ZG + percent_gdp
                             + fuel_income_dependence + VAT.Rate + Crisis.Year
                             + as.factor(country2)
                             + as.factor(country2)*time_trend_country, 
                             data=censdat[censdat$democracy==1,],
                             clusters = country2)

# Summaries for Table 1
texreg(list(fit.lfe1.dem, fit.cfe1.dem, fit.lead.traits.dem, fit.crisis1.dem), 
       digits = 4, omit.coef = c('as.factor(obsid)*|as.factor(country2)*'), include.ci = FALSE)

####################################################
######## Supporting Information -- Table S2 ########                     
####################################################

#Table S2. Model results for non-democracies

#For all models, column 1 has country intercepts and country trends,
#column 2 has country intercepts and leader traits
#column 3 has leader fixed effects

#Table 3: All Autocracies

fit.cfe1.aut <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(country2)
                          + as.factor(country2)*time_trend_country, 
                          data=censdat[censdat$democracy==0,],
                          clusters = country2)

fit.lead.traits.aut <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                                 + NY.GDP.MKTP.KD.ZG + percent_gdp
                                 + fuel_income_dependence + VAT.Rate
                                 + leader_age + gender + unieducation + execrlc
                                 + as.factor(country2)
                                 + as.factor(country2)*time_trend_country,
                                 data=censdat[censdat$democracy==0,],
                                 clusters = country2)

fit.lfe1.aut <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(obsid)
                          + as.factor(obsid)*time_trend_leader, 
                          data=censdat[censdat$democracy==0,],
                          clusters = country2)

fit.crisis1.aut <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                             + NY.GDP.MKTP.KD.ZG + percent_gdp
                             + fuel_income_dependence + VAT.Rate + Crisis.Year
                             + as.factor(country2)
                             + as.factor(country2)*time_trend_country, 
                             data=censdat[censdat$democracy==0,],
                             clusters = country2)

# Summaries for Table 1
texreg(list(fit.lfe1.aut, fit.cfe1.aut, fit.lead.traits.aut, fit.crisis1.aut), 
       digits = 4, omit.coef = c('as.factor(obsid)*|as.factor(country2)*'), include.ci = FALSE)

####################################################
######## Supporting Information -- Table S3 ########                     
####################################################

#Table S3. Model results for parliamentary democracies

#For all models, column 1 has country intercepts and country trends,
#column 2 has country intercepts and leader traits
#column 3 has leader fixed effects

#Table 4: Parliamentary Democracies

fit.cfe1.par <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(country2)
                          + as.factor(country2)*time_trend_country, 
                          data=censdat[censdat$democracy==1 & censdat$parliamentary==2,],
                          clusters = country2)

fit.lead.traits.par <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                                 + NY.GDP.MKTP.KD.ZG + percent_gdp
                                 + fuel_income_dependence + VAT.Rate
                                 + leader_age + gender + unieducation + execrlc
                                 + as.factor(country2)
                                 + as.factor(country2)*time_trend_country,
                                 data=censdat[censdat$democracy==1 & censdat$parliamentary==2,],
                                 clusters = country2)

fit.lfe1.par <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(obsid)
                          + as.factor(obsid)*time_trend_leader, 
                          data=censdat[censdat$democracy==1 & censdat$parliamentary==2,],
                          clusters = country2)

fit.crisis1.par <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                             + NY.GDP.MKTP.KD.ZG + percent_gdp
                             + fuel_income_dependence + VAT.Rate + Crisis.Year
                             + as.factor(country2)
                             + as.factor(country2)*time_trend_country, 
                             data=censdat[censdat$democracy==1 & censdat$parliamentary==2,],
                             clusters = country2)

# Summaries for Table 4
texreg(list(fit.lfe1.par, fit.cfe1.par, fit.lead.traits.par, fit.crisis1.par), 
       digits = 4, omit.coef = c('as.factor(obsid)*|as.factor(country2)*'), include.ci = FALSE)

####################################################
######## Supporting Information -- Table S4 ########                     
####################################################

#Table S4. Model results for presidential democracies

#For all models, column 1 has country intercepts and country trends,
#column 2 has country intercepts and leader traits
#column 3 has leader fixed effects

#Table 5: Presidential Democracies

fit.cfe1.prs <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(country2)
                          + as.factor(country2)*time_trend_country, 
                          data=censdat[censdat$democracy==1 & censdat$parliamentary==0,],
                          clusters = country2)

fit.lead.traits.prs <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                                 + NY.GDP.MKTP.KD.ZG + percent_gdp
                                 + fuel_income_dependence + VAT.Rate
                                 + leader_age + gender + unieducation + execrlc
                                 + as.factor(country2)
                                 + as.factor(country2)*time_trend_country,
                                 data=censdat[censdat$democracy==1 & censdat$parliamentary==0,],
                                 clusters = country2)

fit.lfe1.prs <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(obsid)
                          + as.factor(obsid)*time_trend_leader, 
                          data=censdat[censdat$democracy==1 & censdat$parliamentary==0,],
                          clusters = country2)

fit.crisis1.prs <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                             + NY.GDP.MKTP.KD.ZG + percent_gdp
                             + fuel_income_dependence + VAT.Rate + Crisis.Year
                             + as.factor(country2)
                             + as.factor(country2)*time_trend_country, 
                             data=censdat[censdat$democracy==1 & censdat$parliamentary==0,],
                             clusters = country2)

# Summaries for Table 5
texreg(list(fit.lfe1.prs, fit.cfe1.prs, fit.lead.traits.prs, fit.crisis1.prs), 
       digits = 4, omit.coef = c('as.factor(obsid)*|as.factor(country2)*'), include.ci = FALSE)


####################################################
######## Supporting Information -- Table S5 ########                     
####################################################

#Table S5. Model results for oil importers

#For all models, column 1 has country intercepts and country trends,
#column 2 has country intercepts and leader traits
#column 3 has leader fixed effects

#Table 6: Oil Importers

fit.cfe1.imp <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(country2)
                          + as.factor(country2)*time_trend_country, 
                          data=censdat[censdat$oil==0,],
                          clusters = country2)

fit.lead.traits.imp <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                                 + NY.GDP.MKTP.KD.ZG + percent_gdp
                                 + fuel_income_dependence + VAT.Rate
                                 + leader_age + gender + unieducation + execrlc
                                 + as.factor(country2)
                                 + as.factor(country2)*time_trend_country,
                                 data=censdat[censdat$oil==0,],
                                 clusters = country2)

fit.lfe1.imp <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(obsid)
                          + as.factor(obsid)*time_trend_leader, 
                          data=censdat[censdat$oil==0,],
                          clusters = country2)

fit.crisis1.imp <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                             + NY.GDP.MKTP.KD.ZG + percent_gdp
                             + fuel_income_dependence + VAT.Rate + Crisis.Year
                             + as.factor(country2)
                             + as.factor(country2)*time_trend_country, 
                             data=censdat[censdat$oil==0,],
                             clusters = country2)

# Summaries for Table 6
texreg(list(fit.lfe1.imp, fit.cfe1.imp, fit.lead.traits.imp, fit.crisis1.imp), 
       digits = 4, omit.coef = c('as.factor(obsid)*|as.factor(country2)*'), include.ci = FALSE)

####################################################
######## Supporting Information -- Table S6 ########                     
####################################################

#Table S6. Model results for oil exporters

#For all models, column 1 has country intercepts and country trends,
#column 2 has country intercepts and leader traits
#column 3 has leader fixed effects

#Table 7: Oil Exporters

fit.cfe1.exp <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(country2)
                          + as.factor(country2)*time_trend_country, 
                          data=censdat[censdat$oil==1,],
                          clusters = country2)

fit.lead.traits.exp <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                                 + NY.GDP.MKTP.KD.ZG + percent_gdp
                                 + fuel_income_dependence + VAT.Rate
                                 + leader_age + gender + unieducation + execrlc
                                 + as.factor(country2)
                                 + as.factor(country2)*time_trend_country,
                                 data=censdat[censdat$oil==1,],
                                 clusters = country2)

fit.lfe1.exp <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(obsid)
                          + as.factor(obsid)*time_trend_leader, 
                          data=censdat[censdat$oil==1,],
                          clusters = country2)

fit.crisis1.exp <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                             + NY.GDP.MKTP.KD.ZG + percent_gdp
                             + fuel_income_dependence + VAT.Rate + Crisis.Year
                             + as.factor(country2)
                             + as.factor(country2)*time_trend_country, 
                             data=censdat[censdat$oil==1,],
                             clusters = country2)

# Summaries for Table 7
texreg(list(fit.lfe1.exp, fit.cfe1.exp, fit.lead.traits.exp, fit.crisis1.exp), 
       digits = 4, omit.coef = c('as.factor(obsid)*|as.factor(country2)*'), include.ci = FALSE)

####################################################
######## Supporting Information -- Table S7 ########                     
####################################################

#Table S7. Model results for persistent gasoline subsidizers

#For all models, column 1 has country intercepts and country trends,
#column 2 has country intercepts and leader traits
#column 3 has leader fixed effects

#Table 8: Fuel Persistent Subsidizers

fit.cfe1.sub <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(country2)
                          + as.factor(country2)*time_trend_country, 
                          data=censdat[censdat$persistent==1,],
                          clusters = country2)

fit.lead.traits.sub <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                                 + NY.GDP.MKTP.KD.ZG + percent_gdp
                                 + fuel_income_dependence + VAT.Rate
                                 + leader_age + gender + unieducation + execrlc
                                 + as.factor(country2)
                                 + as.factor(country2)*time_trend_country,
                                 data=censdat[censdat$persistent==1,],
                                 clusters = country2)

fit.lfe1.sub <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(obsid)
                          + as.factor(obsid)*time_trend_leader, 
                          data=censdat[censdat$persistent==1,],
                          clusters = country2)

fit.crisis1.sub <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                             + NY.GDP.MKTP.KD.ZG + percent_gdp
                             + fuel_income_dependence + VAT.Rate + Crisis.Year
                             + as.factor(country2)
                             + as.factor(country2)*time_trend_country, 
                             data=censdat[censdat$persistent==1,],
                             clusters = country2)

# Summaries for Table 8
texreg(list(fit.lfe1.sub, fit.cfe1.sub, fit.lead.traits.sub, fit.crisis1.sub), 
       digits = 4, omit.coef = c('as.factor(obsid)*|as.factor(country2)*'), include.ci = FALSE)

####################################################
######## Supporting Information -- Table S8 ########                     
####################################################

#Table S8. Model results for persistent gasoline taxers

#For all models, column 1 has country intercepts and country trends,
#column 2 has country intercepts and leader traits
#column 3 has leader fixed effects

#Table 9: Fuel Persistent Taxers

fit.cfe1.tax <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(country2)
                          + as.factor(country2)*time_trend_country, 
                          data=censdat[censdat$fixedprice == 1 & censdat$persistent == 0,],
                          clusters = country2)

fit.lead.traits.tax <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                                 + NY.GDP.MKTP.KD.ZG + percent_gdp
                                 + fuel_income_dependence + VAT.Rate
                                 + leader_age + gender + unieducation + execrlc
                                 + as.factor(country2)
                                 + as.factor(country2)*time_trend_country,
                                 data=censdat[censdat$fixedprice == 1 & censdat$persistent == 0,],
                                 clusters = country2)

fit.lfe1.tax <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                          + NY.GDP.MKTP.KD.ZG + percent_gdp
                          + fuel_income_dependence + VAT.Rate
                          + as.factor(obsid)
                          + as.factor(obsid)*time_trend_leader, 
                          data=censdat[censdat$fixedprice == 1 & censdat$persistent == 0,],
                          clusters = country2)

fit.crisis1.tax <- lm_robust(bmgap2015adj ~ loggdppc + I(loggdppc^2)
                             + NY.GDP.MKTP.KD.ZG + percent_gdp
                             + fuel_income_dependence + VAT.Rate + Crisis.Year
                             + as.factor(country2)
                             + as.factor(country2)*time_trend_country, 
                             data=censdat[censdat$fixedprice == 1 & censdat$persistent == 0,],
                             clusters = country2)

# Summaries for Table 9
texreg(list(fit.lfe1.tax, fit.cfe1.tax, fit.lead.traits.tax, fit.crisis1.tax), 
       digits = 4, omit.coef = c('as.factor(obsid)*|as.factor(country2)*'), include.ci = FALSE)


####################################################
######## Supporting Information -- Table S9 ########                     
####################################################

#Table S9. Comparison of rate of change in the net tax across subgroups


mod1 <- lm_robust(rate_change_tax_eff ~ democracy,
                  data=leaders_taxchange_cov,
                  clusters = country2)

mod2 <- lm_robust(rate_change_tax_eff ~ presidential_democracy,
                  data=leaders_taxchange_cov,
                  clusters = country2)

mod3 <- lm_robust(rate_change_tax_eff ~ oil_exporter,
                  data=leaders_taxchange_cov,
                  clusters = country2)

mod4 <- lm_robust(rate_change_tax_eff ~ persistent,
                  data=leaders_taxchange_cov,
                  clusters = country2)


# Summaries for Table 1
texreg(list(mod1, mod2, mod3, mod4), digits = 4, include.ci = FALSE)